//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   Cross_Section.hh
 * \author Jeremy Roberts
 * \date   Jul 8, 2011
 * \brief  Cross_Section class definition.
 * \note   Copyright (C) 2011 Jeremy Roberts. 
 */
//---------------------------------------------------------------------------//
// $Rev::                                               $:Rev of last commit
// $Author::                                            $:Author of last commit
// $Date::                                              $:Date of last commit
//---------------------------------------------------------------------------//


#ifndef CROSS_SECTION_HH_
#define CROSS_SECTION_HH_

#include <vector>
#include <string>

#include "../utilities/DBC.hh"
#include "../utilities/Vector_Lite.hh"

namespace slabtran
{

//===========================================================================//
/*!
 * \class Cross_Section
 * \brief Defines the cross-sections for one material.
 *
 * This class contains the cross-sections for one material.  The full
 * matrix of data is contained in the \ref Materials class.
 *
 * Currently, the total, total absorption, fission, and scattering
 * cross-sections are contained, as are \f$\nu\f$ and \f$\chi\f$.  No
 * (n,g) framework is included, though pure gamma problems could probably
 * be modeled.
 *
 * To represent the scattering cross-section across all groups, it is
 * convenient to use operator notation.  The complete operator is denoted
 * \f$\mathbf{S}\f$, which is defined
 * \f[
 \mathbf{S} = \left(
 \begin{array}{llll}
 \left[\mathbf{S}\right]_{00} & \left[\mathbf{S}\right]_{01}
 & \ldots & \left[\mathbf{S}\right]_{0G} \\
   \left[\mathbf{S}\right]_{10} & \left[\mathbf{S}\right]_{11}
 & \ldots & \left[\mathbf{S}\right]_{1G} \\
   \vdots & \ddots & \ddots & \vdots \\
   \left[\mathbf{S}\right]_{G0} & \left[\mathbf{S}\right]_{G1}
 & \ldots & \left[\mathbf{S}\right]_{GG} \\
   \end{array}\right)\:,
 * \f]
 * where the subscripts imply \f$g\gets g'\f$.  One block
 * \f$[\mathbf{S}]_{gg'}\f$ is a matrix consisting of the single
 * (major) diagonal
 * \f[
 \Big \{ \Sigma_s[l_0,g\gets g'], \, \Sigma_s[l_1,g\gets g'], \,
 \Sigma_s[l_1,g\gets g'], \, \Sigma_s[l_1,g\gets g'], \,
 \Sigma_s[l_1,g\gets g'], \, \Sigma_s[l_2,g\gets g'], \, \ldots , \,
 \Sigma_s[l_L,g\gets g'] \Big \} \, ,
 * \f]
 * given here for 3d, where we note a given Legendre order is repeated for each spherical
 * harmonic with the given \f$l\f$.  These repeated values are not stored.
 *
 * A single row of
 * \f$\mathbf{S}\f$ is used to define the single group problem
 * \f[
 \mathbf{L}[\psi]_g = [\mathbf{M}]([\mathbf{S}]_{g0}[\phi]_0
 + [\mathbf{S}]_{g1}[\phi]_1 + [\mathbf{S}]_{g2}[\phi]_2 +
 \ldots + [\mathbf{S}]_{gG}[\phi]_G) + [\mathbf{M}][q_e]_g \, ,
 * \f]
 * where \f$[\mathbf{\phi}]_g \f$ are the flux moments for a cell for a group.
 *
 * In general, \f$\mathbf{S}\f$ is not a full matrix.  To avoid storing
 * zeros and to facilitate generation of the scattering sources, we store
 * only nonzero values and define for each row and column the lower and
 * upper group indices that contain the nonzero data.  The column bounds
 * correspond to inscatter and are used in typical forward calculations.
 * The row bounds correspond to outscatter, needed for adjoint calculations.
 *
 * In other words, for each row, we store the lower
 * and upper group bounds for in-scatter, i.e. the lowest group \f$g'\f$ that
 * yields downscatter to \f$g\f$ and the highest group \f$g'\f$ that yields
 * upscatter into \f$g\f$.  For each column \f$g\f$, we store the lowest group
 * \f$g'\f$ to which \f$g\f$ can upscatter, etc.
 *
 * All cross sections are in units of \f$\mbox{cm}^{-1}\f$.
 *
 * \todo is it worth computing/storing bounds for fission groups, i.e. chi>0?
 * \todo checks for correct xs, like sum of all = total, sum(chi)=1, etc.
 *
 */
/*! \example test/testCross_Section.cc
 *
 * Test of Cross_Section.
 */
//===========================================================================//

class Cross_Section
{
public:
  /// \name Typedefs
  //\{
  typedef std::vector<double> S_Block;
  typedef std::vector<S_Block> S_Row;
  typedef std::vector<S_Row> S_Operator;
  typedef std::vector<int> Vec_Int;
  typedef std::vector<double> Group_Data;
  //\}

  // Constructor.
  Cross_Section(const int number_groups, const int legendre_order);

  // >>> ACCESSORS/MODIFIERS

  /// Return number of groups.
  int number_groups() const { return d_number_groups; }

  /// Return Legendre order.
  int legendre_order() const { return d_legendre_order; }

  /// Return downscatter flag.
  bool downscatter() const { return d_downscatter; }

  /// Return adjoint.
  int adjoint() const { return d_adjoint; }

  /// Set adjoint
  void set_adjoint(int adjoint) { d_adjoint = adjoint; }

  /// Set name
  void set_name(std::string name) { d_name = name; }

  /// Get name
  std::string name() const { return d_name; }

  //-------------------------------------------------------------------------//

  /// Mutable access of total group cross section in /cm.
  double& sigma(int g) { return d_sigma[g]; }

  /// Const access of total group cross section in /cm.
  double sigma(int g) const { return d_sigma[g]; }

  /// Mutable access of total absorption group cross section in /cm.
  double& sigma_a(int g) { return d_sigma_a[g]; }

  /// Const access of total absorption group cross section in /cm.
  double sigma_a(int g) const { return d_sigma_a[g]; }

  /// Mutable access of total absorption group cross section in /cm.
  double& sigma_f(int g) { return d_sigma_f[g]; }

  /// Const access of fission group cross section in /cm.
  double sigma_f(int g) const { return d_sigma_f[g]; }

  /// Mutable access of nu group value
  double& nu(int g) { return d_nu[g]; }

  /// Const access of nu group value.
  double nu(int g) const { return d_nu[g]; }

  /// Mutable access of chi group value.
  double& chi(int g) { return d_chi[g]; }

  /// Const access of chi group value.
  double chi(int g) const { return d_chi[g]; }

  /*
   *  \brief Mutable access of scattering cross section moment in /cm.
   *
   *  For the forward problem
   *  \f[
   *      \mathbf{T}\psi_g = \mathbf{M}\mathbf{S}_{g\gets 1} +
   *          \mathbf{M}\mathbf{S}_{g\gets 2} + \ldots + S_g \f$
   *  \f]
   *  we need to return a $\Sigma_{g\gets g'}$.  For the adjoint
   *  problem
   *  \f[
   *      \mathbf{T}^*\psi^*_g = \mathbf{M}\mathbf{S}_{g\to 1} +
   *          \mathbf{M}\mathbf{S}_{g\to 2} + \ldots + S^*_g \f$
   *  \f]
   *  we need to return a $\Sigma_{g\to g'}$.  However, this is
   *  just $\Sigma_{g' \from g}$, i.e. we need to switch the
   *  given group indices to access the right data.
   *
   *  \param g        group associated with group equation
   *  \param g_prime  group from(to) which I receive(scatter)
   */
  double& sigma_s(int g, int g_prime, int l);

  /*!
   *  \brief Const access of scattering cross section moment in /cm.
   *
   *  \param g        group associated with group equation
   *  \param g_prime  group from(to) which I receive(scatter)
   */
  const double& sigma_s(int g, int g_prime, int l) const;

  /*!
   *  \brief Mutable access to block of \f$\mathbf{S}\f$.
   *
   *  \param g        group of group equation
   *  \param g_prime  group from(to) receive(scatter)
   */
  std::vector<double>& sigma_s(int g, int g_prime);

  /*!
   *  \brief Const access to block of \f$\mathbf{S}\f$.
   *
   *  \param g        group of group equation
   *  \param g_prime  group from(to) receive(scatter)
   */
  const std::vector<double>& sigma_s(int g, int g_prime) const;

  /// Mutable access to lower bound for scattering loop
  int& group_lower_bound(int g, int adj){return d_group_lower_bound[adj][g];}

  /// Mutable access to upper bound for scattering loop
  int& group_upper_bound(int g, int adj){return d_group_upper_bound[adj][g];}

  /// Lower bound for scattering loop
  int group_lower_bound(int g) const;

  /// Upper bound for scattering loop
  int group_upper_bound(int g) const;

  /// Print the cross-section
  void print() const;

  //-------------------------------------------------------------------------//

private:

  /// Number of groups.
  int d_number_groups;

  /// Legendre order of scattering.
  int d_legendre_order;

  /// Boolean flag for downscatter-only.
  bool d_downscatter;

  /// Flag for adjoint
  int d_adjoint;

  /// Material description
  std::string d_name;

  //-----------------------------------------------------------------------//
  /*!
   * \name Group Data
   *
   * Add additional parameters if needed (e.g (n,2n), gamma's...)
   *
   */
  //\{

  /// Total cross section.
  Group_Data d_sigma;

  /// Absorption (includes radiative capture, fission, etc.) cross section
  Group_Data d_sigma_a;

  /// Fission cross section
  Group_Data d_sigma_f;

  /// Fission neutron yield
  Group_Data d_nu;

  /// Fission (\f$\chi\f$) spectrum
  Group_Data d_chi;

  /// Group-to-group blocks of scattering moments.
  S_Operator d_S;

  /// Group bounds (eV)
  // \todo add me
  std::vector<double> d_group_bounds;

  //\}
  //-----------------------------------------------------------------------//

  /*!
   *  \brief Lower bound for scattering.
   *
   *  Forward: Lowest g that downscatters to me.
   *  Adjoint: Lowest g to which I upscatter.
   */
  util::Vector_Lite<Vec_Int,2> d_group_lower_bound;

  /*!
   *  \brief Upper bound for scattering.
   *
   *  Forward: Highest g that upscatters to me.
   *  Adjoint: Highest g to which I down.
   */
  util::Vector_Lite<Vec_Int,2> d_group_upper_bound;

  // >>> IMPLEMENTATION

  /// Get correct indices for adjoint
  void group_index(int &g1, int &g2) const;

};

} // end namespace slabtran

//---------------------------------------------------------------------------//
// INLINE FUNCTIONS
//---------------------------------------------------------------------------//

#include "material/Cross_Section.i.hh"

#endif /* CROSS_SECTION_HH_ */

//---------------------------------------------------------------------------//
//              end of Cross_Section.hh
//---------------------------------------------------------------------------//
